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Abstract —This paper considers future distribution networks 
featuring inverter-interfaced photovoltaic (PV) systems, and ad¬ 
dresses the synthesis of feedback controllers that seek real- and 
reactive-power inverter setpoints corresponding to AC optimal 
power flow (OPF) solutions. The objective is to bridge the 
temporal gap between long-term system optimization and real¬ 
time inverter control, and enable seamless PV-owner participa¬ 
tion without compromising system efficiency and stability. The 
design of the controllers is grounded on a dual e-subgradient 
method, and semidefinite programming relaxations are advocated 
to bypass the non-convexity of AC OPF formulations. Global 
convergence of inverter output powers is analytically established 
for diminishing stepsize rules for cases where: i) computational 
limits dictate asynchronous updates of the controller signals, and 
ii) inverter reference inputs may be updated at a faster rate than 
the power-output settling time. 

Index Terms —Distribution systems, photovoltaic inverter con¬ 
trol, distributed optimization and control; optimal power flow. 


I. Introduction 

Present-generation residential photovoltaic (PV) inverters 
typically operate in a distributed and uncoordinated fashion, 
with the primary objective of maximizing the power extracted 
from PV arrays. With the increased deployment of behind- 
the-meter PV systems, an upgrade of medium- and low- 
voltage distribution-system operations and controls is required 
to address emerging efficiency, reliability, and power-quality 
concerns ffl, CD. To this end, several architectural frameworks 
have been proposed for PV-dominant distribution systems 
to broaden the objectives of inverter real-time control, and 
enable inverters to partake in distribution-network optimization 
tasks ||3l-ll6l. 

Past works have addressed the design of distributed real¬ 
time inverter-control strategies to regulate the delivery of real 
and reactive power based on local measurements, so that 
terminal voltages are within acceptable levels 0, 0. On a 
different time scale, centralized and distributed optimal power 
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flow (OPF) algorithms have been proposed to compute optimal 
steady-state inverter setpoints, so that power losses and voltage 
deviations are minimized and economic benefits to end-users 
providing ancillary services are maximized 0, 0-EU. 

In an effort to bridge the temporal gap between real¬ 
time control and network-wide steady-state optimization, this 
paper addresses the synthesis of feedback controllers that seek 
optimal PV-inverter power setpoints corresponding to AC OPF 
solutions. The guiding motivation is to ensure that PV-system 
operation and control strategies are adaptable to changing 
ambient conditions and loads, and enable seamless end-user 
participation without compromising system efficiency. The 
proposed feedback controllers continuously pursue solutions 
of the formulated OPF problem by dynamically updating 
the setpoints based on current system outputs and problem 
parameters. This presents significant improvements over state- 
of-the-art distributed OPF approaches wherein reference sig¬ 
nals for the inverters are updated only upon convergence of 
the distributed algorithm. In this setting, it is evident that if 
problem parameters or inputs change during the computation, 
broadcast, and implementation of the distributed solution of 
the OPF, the inverter would implement outdated setpoints. 

Prior efforts in this direction include continuous-time feed¬ 
back controllers that seek Karush-Kuhn-Tucker conditions for 
optimality developed in fl2l . and applied to solve an economic 
dispatch problem for bulk power systems in E|. Recently, 
modified automatic generation and frequency control methods 
that incorporate optimization objectives corresponding to DC 
OPF problems have been proposed for lossless bulk power 
systems in e.g., HD-Ill). A heuristic based on saddle-point- 
flow methods is utilized in El to synthesize controllers 
seeking AC OPF solutions. Strategies that integrate economic 
optimization within droop control for islanded lossless micro¬ 
grids are developed in ffl8l . In a nutshell, these approaches 
are close in spirit to the seminal work ED, where dynamical 
systems that serve as proxies for optimization variables and 
multipliers are synthesized to evolve in a continuous-time 
gradient-like fashion to the saddle points of the Lagrangian 
function associated with a convex optimization problem. For 
DC OPF, a heuristic comprising continuous-time dual ascent 
and discrete-time reference-signal updates is proposed in (20j; 
where, local stability of the resultant closed-loop system is also 
established. 

Distinct from past efforts jl3l - llT8l . Il20l , this work lever¬ 
ages dual e-subgradient methods ED, E2, to develop a 
feedback controller that steers the inverter output powers 
towards the solution of an AC OPF problem. A semidefinite 
programming (SDP) relaxation is advocated to bypass the non- 
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convexity of the formulated AC OPF problem IfTOl , 1231 , 124S . 
The proposed scheme involves the update of dual and primal 
variables in a discrete-time fashion, with the latter constituting 
the reference-input signals for the PV inverters. Convergence 
of PV-inverter-output powers to the solution of the formulated 
OPF problem is analytically established for settings where: 
i) in an effort to bridge the time-scale separation between 
optimization and control, the reference inputs may be updated 
at a faster rate than the power-output settling time; and, ii) 
due to inherent computational limits related to existing SDP 
solvers, the controller signals are updated asynchronously. 
Although the present paper focuses on the case where an 
SDP relaxation is utilized to bypass the non-convexity of 
OPF tasks, the proposed synthesis procedure can be utilized 
to develop controllers that provably drive the inverter output 
to solutions of various convex relaxations |251 and linear 
approximations I26l - lf28fl of the OPF problem. 

Overall, the proposed framework considerably broadens the 
approaches of fT3l - fT8l . {20) by: i) considering AC OPF 
setups; ii) incorporating PV-inverter operational constraints; 
Hi) accounting for communication constraints which naturally 
lead to discrete-time controller updates; and, vi) accounting for 
computational limits which involves an asynchronous update 
of the controller signals. It is also shown that the controller 
affords a distributed implementation, and requires limited 
message passing between the PV systems and the utility. 

The remainder of this paper is organized as follows. 
Section [IT] outlines the problem formulation, while the PV 
controller is developed in Section [HI] Section [IV] elaborates 
on the distributed implementation of the proposed control 
architecture. Numerical tests are reported in Section [V] and 
conclusions are provided in Section (VT] 


II. Problem Formulation 


Dynamical models and relevant formulations for optimiz¬ 
ing inverter setpoints are outlined for a general networked 
dynamical system in Section III-A 


PV-inverter control in Section III-B 


and tailored to real-time 


1 Notation. Upper-case (lower-case) boldface letters will be used for ma¬ 
trices (column vectors); (-) T for transposition; (•)* complex-conjugate; and, 
(•) H complex-conjugate transposition; 5R{-} and ■ \ denote the real and 
imaginary parts of a complex number, respectively; j := y/—l. Tr(-) the 
matrix trace; rank(-) the matrix rank; | • | denotes the magnitude of a number or 
the cardinality of a set; vec(X) returns a vector stacking the columns of matrix 
X, and bdiag({Xi}) forms a block-diagonal matrix. R ;V and C N denote 
the spaces of N X 1 real-valued and complex-valued vectors, respectively; N 
the set of natural numbers; and, W+ xN denotes the space of N X IV positive 
semidefinite Hermitian matrices. Given vector x and square matrix X, | |x ||2 
denotes the Euclidean norm of x, and || V|| 2 the (induced) spectral norm of 
matrix X. [x]j ([(x)]i) points to the j-th element of a vector x (vector-valued 
function f(x)). x(t) is the time derivative of x(t). Given a scalar function 
/(x) : R™ —> R, V x /(x) returns the gradient [J ^,..., ^-] T . For a 
continuous function /(£), /[£*.] denotes its value sampled at t Finally, Ijy 
denotes the N xN identity matrix; and, OmxN> ImxJV the MxN matrices 
with all zeroes and ones, respectively. 


A. General problem setup 

Consider Nd dynamical systems described by 

±i(t) = £i (x,;(f), dj(f), Ui (f)) (la) 

Yi{t) = r i (x i (t),d i (t)'), i G N d ■= {1,...,N D } (lb) 

where: Xj(i) G R"' 1 * is the state of the i-th dynamical 
system at time t; yj(f) £ ^ C R"^ is the measurement 
of state Xj(f) at time t\ u,(f) G Vi is the reference input; 
and dj(i) G Vi C R"'*’ 4 is the exogenous input. Finally, 
f i : R ” x ’ 4 x R”*R x R” tl . i x R"vR —► R” x R and ir : 
r«*,; x _y R n »,< are arbitrary (non)linear functions. The 

following system behavior for given finite exogenous inputs 
and reference signals is assumed. 

Assumption 1: For given constant exogenous inputs {d, G 
^iliGA/ri an d reference signals {Uj G Vi}iejV D , there exist 
equilibrium points {xi}i e yv D for ([]} that satisfy: 

0 = f» (xj, dj, Uj) (2a) 

Ui = r i (x i; dj) , i G Jf D ■ (2b) 

Notice that in (l2bl > the equilibrium output coincides with the 
commanded input u,; that is, y,; = u,. These equilibrium 
points are locally asymptotically stable ll29l . □ 

For given exogenous inputs {d, G 'Z2,;}, p _,v'd, consider the 
following optimization problem: 

(PI) min H{V) + Y Gi(ui) (3a) 

vev.iu.ey,} 

subject to hj(V) — Uj + dj = 0, V* G Mo (3b) 

where V C H"' /Xrtv is a convex, closed, and bounded subset 
of the cone of positive semidefinite (Hermitian) matrices; 
function H(V) : i”' ,xnv —> R is known, convex and finite 
over V; Gj(u,) : R rl “- i — > R is strongly convex and finite over 
Vi; and, the vector-valued function h, ; (V) : M((. vXnv —> R"v,< 
is affine. Finally, sets {Vi}iGA/" D > which define the space 
of possible reference inputs for the dynamical systems, are 
assumed to comply to the following requirement. 

Assumption 2: Sets {Vi}ieA/b are convex and compact. 
Further, (PI) has a non-empty feasible set and a finite optimal 
cost. □ 

With these assumptions, problem (PI) is a convex program; 
moreover, it can be reformulated into a standard SDP form by 
resorting to the epigraph form of the cost function. 

It is evident from (l2bl i that (PI) defines the optimal op¬ 
erating setpoints of the dynamical systems CD in terms of 
steady-state outputs CO, ESI- In fact, by utilizing the optimal 
solution {u° pt }ig7v' D °f (PI) as reference inputs, it follows 
from (l2bl > that each system output will eventually be driven to 
the point y, : = u° pt . Function (l3at captures costs incurred by 
the steady-state outputs, as well as costs associated with matrix 
variable V, which couples the steady-state system outputs 
{y i = u i}iej\r D through the linear equality constraints (l3bl >. 

In principle, (PI) could be solved centrally by a system- 
level control unit, which subsequently dispatches the reference 
signals {u° pt }i ej y D for the dynamical systems. In lieu of a 
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centralized solution of (PI), the objective here is to design a 
distributed feedback controller for the dynamical systems ©. 
so that the resultant closed-loop system is globally convergent 
to an equilibrium point {x,;} ieJ v D , {y t = r;(x;, dj)} ieJ y D , 
where the values {yi}i£j\f D of the steady-state outputs coin¬ 
cide with the optimal solution {u° pt } ieNo of ( P1 )- 


B. PV-inverter output regulation to OPF solutions 

The task of regulating the power output of PV inverters is 
outlined in this section, and cast within the framework of 0 - 
0. In this regard, 0-0 will model the inverter dynamics lf30l 
Ch. 8], PTI : while OPF will be formulated in the form 0 by 
leveraging SDP relaxation techniques l23l . l24l . 

Network. Consider a distribution system comprising TV + 1 
nodes collected in the set A f, and lines represented by the set 
of undirected edges £ := {(to, n) : to, n £ M}. The set TV" := 
{0,1,..., TV} is partitioned as TV" = {0} U Md U TV"o, where: 
node 0 denotes the secondary of the step-down transformer; 
inverter-interfaced PV systems are located at nodes TVd = 
{1,..., TVd} [cf. 0 ]; and, Mo '■= {Md + 1,..., TV} collects 
nodes with no power generation. For simplicity of exposition, 
the framework is outlined for a balanced system; however, the 
proposed framework can be extended to unbalanced multi¬ 
phase systems as explained in Appendix iDl 


Dynamics of PV inverters. Equation ( flat is utilized to model 
the dynamics of PV inverters, regulating real- and reactive 
output powers to prescribed setpoints. For example, relevant 
dynamical models for inverters operating in a grid-connected 
mode are discussed in e.g., If30l Ch. 8] and PTi . These models 
can be conveniently cast within 0-0 as shown next. 

• Let pft) := Eft)cos(cot + fft))ift) and qft ) := 
Eft) cos(yjt +fft) — it/2)ift) denote the instantaneous out¬ 
put real and reactive powers of inverter i £ Md, respectively, 
where w is the grid frequency, vft) := Eft) cos (cot + 4>ft)) 
the voltage waveform, and ii (t) is the current injected. Further, 
let Pi (t) and Qi (t) denote averages of the instantaneous output 
real and reactive powers over an AC cycle; that is, 


Pi(t) ■= 


UJ 
27r 



Qft) 


UJ 

2tt 



(4) 


Then, the state of system 0 is ~x.ft) := [Pft), Qft)f. 

• Vector u, (t) = u, collects the constant commanded real and 
reactive powers for inverter i. By 0 , inverters regulate the 
output powers to the commanded setpoints u, ; see e.g., l30l 
Ch. 8], E). 

• Let Pe,ft) and Qeft) denote the demanded real and reactive 
loads at node i £ M. Then, vector d, (t) is set to be dj(f) := 
[Pe,i(t),Qc,ft)] J for all i £ TV\{0}. 

• By setting 

r fxft),dft)) = xft) (5) 


< 1 1 b| i equates the state with the measurement of the inverter 
output powers. 

Steady-state OPF problem. Let V) := (£ , j/v / 2)e^ i £ C be 
the phasor representation of the steady-state voltage at node 
i £ TV". Similarly, let £ C denote the phasor for the current 
injected at node i £ M, and define i := [Iq, . . . , In] J £ C Ar+1 



Fig. 1. Operating regions 3^ for PV inverter under: (a) reactive power 
compensation 0; (b) real power curtailment 0; and, (c) combined real- and 
reactive-power control 0. 


and v := [Vo,..., V;v] T £ C N+1 . Then, using Ohm’s and 
Kirchhoff’s circuit laws, the linear relationship i = Yv can be 
established, where Y £ C N+lxN+1 is the admittance matrix 
formed based on the distribution-network topology and the n- 
equivalent circuits of lines £ ll32l Ch. 6]; see e.g., iflOl . l23l . 
Il24l . |33l for details on the construction of matrix Y. 

For prevailing ambient conditions, let Pf > 0 denote the 
available real power for the inverter i £ Md- The available 
power is a function of the incident irradiance, and corresponds 
to the maximum power point of the PV array. When PV- 
systems operate at unity power factor El, a set of challenges 
related to power quality and reliability in distribution sys¬ 
tems may emerge for sufficiently high levels if deployed PV 
capacity 0 . For instance, overvoltages may be experienced 
during periods when PV generation exceeds the household 
demand 0 , while fast-variations in the PV-output tend to 
cause transients that lead to wear-out of legacy switchgear 121. 
Efforts to ensure reliable operation of existing distribution 
systems with increased behind-the-meter PV generation are 
focused on the possibility of inverters providing reactive power 
compensation and/or curtailing real power 0-El. In the most 
general setting, the set of operating points for PV inverters 
providing ancillary services can be specified as: 


Vi = {(Pi,Qi)'- Pr <Pi< p?, Qi < s 2 i - pf, 

and |Qi| < (tan 9)Pi] (6) 


where Si is the rated apparent power, and the last inequality 
is utilized to enforce a minimum power factor of cos 9. 
Parameters 9 and P™ n can be conveniently tuned to account 
for the following strategies: 

(cl) Reactive power compensation: P™ n = Pf, 6 £ (0,7r/2]; 

(c2) Active power curtailment: p™ n £ [0,P? V ), 9 = 0; and, 

(c3) Active and reactive control: P™ n £ [0, P? v ], 9 £ (0, 7t/2]. 

The PV-inverter operating regions involved by strategies (cl)- 
(c3) are illustrated in Figure 0 It is evident that sets {Vi} 
adhere to Assumption [2] 

For given load and ambient conditions, a prototypical OPF 
formulation for optimizing the steady-state operation of a 
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distribution system is given as follows: 

(OPF) min H(v)+ Y G l (P ll Q l ) (7a) 

vMP " Qi} i£* D 

subject to i = Yv, and 

Vil* = P t - P tli + j (Qi - Q t ,i), Vi e M d (7b) 
V n I* = -P tin - j Q t , n , VneMo (7c) 

V min < \Vi\ < V max Vi g M (7d) 

Uieji; Vie A/d (7e) 


where l /,nln and l /,n:ix are prescribed voltage limits (e.g., 
ANSI C84.1 limits); the constraint on | Vo | is left implicit; (ITel i 
specifies the feasible inverter operating region [cf. Figure [Tj; 
and, equalities (ITbll -dTcl) capture the power-balance equations 
for nodes with and without inverters, respectively. For nodes 
without loads (e.g., utility poles), one clearly has that P(- i = 

Ql,i = 0. 

Function H (v) can capture various network-oriented perfor¬ 
mance objectives that the distribution system operator (DSO) 
may pursue. For example, the DSO may aim to minimize 
the power losses on the distribution lines, voltage magnitude 
deviations from nominal, and/or the power drawn from the 
substation 0, ED- On the other hand, function Gi(Pi,Qi ) 
models PV-inverter costs/rewards for ancillary service provi¬ 
sioning such as real power curtailment and/or reactive power 
compensation 0, 0. ED; for example, this function can be 
set to Gi(Pi,Qi) = afPr - F) 2 + h{P™ - Pi) + CiQ? + 
di\Qi\, with di,bi,Ci, di denoting market-oriented coefficients, 
to maximize the amount of power provided by PV systems. 
Finally, notice that additional constraints such as thermal limits 
may be naturally accommodated in ([TJ. 

It is well-known that the OPF problem 0 is nonconvex, 
and thus hard to solve to global optimality in both centralized 
and distributed setups. Further, given that the problem is 
nonconvex, convergence of distributed algorithms (derived, 
e.g., via Lagrangian decomposition techniques) is not always 
guaranteed and needs to be established. Since the objective of 
this work is to develop distributed controllers so that inverter 
output powers are pmvably convergent to OPF solutions, a 
convex reformulation of the OPF task is considered next. 

SDP relaxation of the OPF problem. To formulate an 
SDP relaxation of the prototypical steady-state OPF prob¬ 
lem 0, consider expressing steady-state powers and voltage 
magnitudes as linear functions of the outer-product matrix 
V := vv H (23l , Il24l . and define matrix Y, := e^ejY per 
node i, where denotes the canonical basis of R^M 

Using Y,;, form the Hermitian matrices dy := 4(Y, + Y^), 
:= i(Yj — Y ? H ), and Y, := e,ej. Then, the balance 
equations for real and reactive powers at node i £ Ay.) can be 
expressed as Tr(«l?jV) = Pj — P( t and Tr(\t,V) = Qi —Qn j, 
respectively. To reformulate the OPF in the form 0, consider 
setting Uj = [P i ,Q l ] J , d* = [P<?,j, Qi] T , and 


( 8 ) 


Additionally, define the following convex set: 

V := {V : V h 0, V* in < Tr(YjV) < V^Vi £ M 
and Tr($ 4 V) = -P t% , Tr('iyV) = -Q u Vi e Mo} ■ (9) 

With these definitions, problem 0 can be equivalently ex¬ 
pressed as follows 

min iT(V)+ Gi(uj) (10a) 

vev , {u , ei ) i} ^ 

subject to hj(V) — + d; = 0, Vi e Md (10b) 

rank(V) = 1. (10c) 

On par with 0 , problem (fTOt is nonconvex because of the 
rank constraint; however in the spirit of the SDP relaxation, 
the constraint m can be dropped. Notice that, once the 
constraint (1 1 Ocb is dropped, the resultant SDP relaxation of 
the OPF problem is in the form 0 . If the optimal matrix 
V opt of the relaxed problem has rank(V opt ) = 1, then the 
resultant power flows are globally optimal ED, EM- Sufficient 
conditions for this relaxation to be exact for radial and 
balanced systems are provided in ll34l . while its applicability 
to unbalanced multiphase systems is investigated in Col. 

In this setup, the objective of the feedback controller that 
will be designed in the following section, is to drive the 
inverter outputs {jy(/) = [Pi{t), Qi(/)] T }ieA/b t0 the optimal 
solution {u° p )ieN' r> of the OPF problem. 

III. Feedback Controller 
Dual e-subgradient methods are leveraged in Section IIH-BI 
to synthesize controllers for systems 0 whose outputs track 
recursive solvers of (PI). Applications to the real-time PV- 
inverter control problem are discussed in Section [IV] 

To streamline proofs of relevant analytical results, it will be 
convenient to express the linear equality constraints (l3bl > in a 
compact form. To this end, define u := [uj,... , d := 

[di,.... dT J T and h(V) := [hi(V),..., (V)] T . Then, 

constraints (TThb can be compactly expressed as h(V) = u— d. 


A. Primer on dual gradient method 

Consider the Lagrangian corresponding to 0, namely: 

L (V, {uj}, {Aj}) := iJ(V)+ ^ Gi{ui) 

+ y AjthiCvj-ui+do (id 

where G denotes the Lagrange multiplier associated 
with (l3bl >. Based on CD, the dual function and the dual 
problem are defined as follows (see, e.g., E3): 




min 

VGV,{u ieyi}i GJ * D 


L(V,{ iiil.IM) 


( 12 ) 


5 opt := max g({A,}). (13) 

{Ai} ie A r D 

Regarding the optimal Lagrange multipliers, the following 
technical requirement is presumed in order to guarantee their 
existence and uniqueness; see e.g., (36). 


MV) = [Tr($,V),Tr(^V)] T . 




IEEE TRANSACTIONS ON POWER SYSTEMS 


5 


Assumption 3: Vectors 

V[ vec T(v),uT]T[h(V) +g(u,d)]i, i = l,..., Y n y,i ( 14 ) 

are linearly independent. □ 

Section m will elaborate on how condition © can be 
checked in the OPF context. Under current modeling as¬ 
sumptions, it follows that the duality gap is zero, and the 
dual function g({A;}) is concave, differentiable, and it has 
a continuous first derivative l37l . Consider then utilizing a 
gradient method to solve the dual problem, which amounts to 
iteratively performing l37l : 

{ V[Al], {u,;[/c]}igA/b} 

= arg min L(V, {uj, {A, [A;]}) (15a) 

vev.luje^i} 

A i[k + 1] = A ilk) + a k+1 V x L(V[k], {u,[fc]}, {AJ) (15b) 

where k £ N denotes the iteration index, a k +i > 0 is the 
stepsize, and ( 1 1 5bb is repeated for all i £ J\f. In particular, 
a non-summable but square-summable stepsize sequence is 
adopted in this paper f 22 l : that is, there exist sequences 
{ 7 fc}fc>o and {r] k } k > 0 such that: 

(si) 7 * —>■ 0 as k —► + 00 , and J2 k =o 7 fc = +°° ; 

(s 2 ) 7 fc < a k < rj k for all k > 0 ; and, 

(s3) r/ k i 0 as k -»• +oo, and J2 k =o < +°°- 
At iteration k, the same step-size a k is utilized for all i £ Af. 
Exploiting the decomposablility of the Lagrangian, steps (fT5l > 
can be equivalently expressed as: 

A i[k + 1] = A i[k] + a k+ i (h,(V[£;]) - u t [k\ + d*) (16a) 

u i[k + 1] = arg min Gi(iij) — Aj[fc + l]u* (16b) 

u 

V[k + 1] = arg min H (V) + Y A I[ fc + X ] M V ) (16c) 

i£A/b 

with (1 1 6bb — (1 1 6ab performed for all i £ ASd . Finally, notice 
that from the compactness of sets V and it follows 

that there exists a scalar G, 0 < G < +oo, such that 

||h(V[fc]) — u[fc] + d )|| 2 < G, VfceN. (17) 

Using ©, and a stepsize sequence {a k } k >o satisfying 
(sl)-(s3), it turns out that the dual iterates A,; [k] converge 
to the optimal solution A° pt of the dual problem (IT3l >: that 
is, || A opt - A[/c] || 2 ^ 0 as k —> oo 051 Prop. 8.2.6], I©, 
B 2 - Iterates V[fc] and {tq[fc]}ieA/'D become asymptotically 
feasible and their optimal values, V opt and {u“ pt }i EJ v D , 
can be recovered from (1 16cb and (116bb . respectively, once 
{A° pt }ieA/D becomes available. 

Steps similar to © are typically adopted to enable a 
distributed solution of the OPF 0, 133), [1381-1411 as well 
as other resource allocation tasks such as the economic dis¬ 
patch problem and residential load control l42l . As illustrated 
in Figure [2a) and explained in detail in Section [IV] up¬ 
dates (1 1 6ab - (l 1 6bb are implemented at each individual PV sys¬ 
tem while ( 1 1 6cb is performed at the DSO. However, in conven¬ 
tional approaches, the optimal reference signals {u ( j !pl },c_,v' D 
are implemented at the PV-inverters only when the distributed 
algorithm converges to the optimal solution. It is evident that 


under this operating paradigm the optimization and local con¬ 
trol tasks operate at two different time scales, with reference 
signals updated every time that the OPF problem is solved and 
implemented only when the inverter dynamics are in steady 
state. This motivates the development of control schemes 
that continuously pursue solutions of the OPF problem by 
dynamically updating the setpoints, based on current system 
outputs and problem parameters. This is described next. 

B. Controller synthesis 

Consider updates performed at discrete time instants t £ 
{t k ,k £ N}, with V[t fe ], {u ?; [f fe ]} ieA /- D , and let {Ai[f fc ]} ieA /- D 
denote the values of primal and dual variables, respectively, 
at time t k . The following method accounts for the system 
dynamics in ([Qi while solving (PI) with dual-gradient-based 
approaches. 

At time t k , the system outputs are sampled as: 

y»[ffc] = rj(xi(i fc ),di) ffi£Af D (18a) 

and, they are utilized to update the dual variables as follows: 

[^/c+l] = 

+ a k +i (hj(V[4]) - y i[t k \ + d,J , Vi £ Af D - (18b) 

Given A,[ffc+i], the primal variables V[£fc+i] and 
{u,;[ffc + i]}ig 7 v' D are then updated as: 

Ui[ffc+t] = arg min Gi(u ; ) - [£ fc+ i]u, ; . (18c) 

u i&y>i 

V[£ fc+1 ] = arg nun 7T(V) + Y A I[^+i] h i( v ) (l gd ) 

Once (fl8ct is solved, the vector-valued reference signal 
Uj[ffc + i] is applied to the dynamical system (flat over the inter¬ 
val {t k ,t k + 1 ]; that is, u i(t) = Ui[t k+1 ],t£ (t k ,t k+ 1 ]. At time 
t k +i the system outputs {yi[ffc+i]}ieA/b are sampled again, 
and (U8bb - (118cb are repeated. Notice that, differently from 
standard dual gradient methods, variable u, [t k ] is replaced by 
the sampled system output y i[t k ] in the ascent step (1 1 8bb . 

Steps (1 1 8bb — (1 1 8cb in effect constitute the controller for 
the dynamical systems 0. Specifically, the (continuous-time) 
reference signals {u,;(f)},; e // D produced by the controller have 
step changes at instants {t k ,k £ N}, are left-continuous 
functions, and take the constant values {uifffc+ijjjgjvy, over 
the time interval (t k , t k + 1 ]. It is evident that if u i[t k \ converges 
to u° pt as k —> oo (and thus u,(f) — y u° pt as t —> oo), then 
y i ( t ) —► u° pt as t —> oo by virtue of Q. 

Suppose for now that the interval (t k -i,t k \ is large enough 
to allow the outputs {y.j(f)}j e _y D to converge to the com¬ 
manded input {u,[7,:| }i6.Vd [cf- ©]■ Under this ideal setup 
with a pronounced and tangible time-scale separation between 
controller and system dynamics, one has that lim t ^ t - ||yi(f) — 
u i[ffc]|| = 0, for all k [cf. <0], and step (1 1 8bb is replaced by 
Ai[tk+i] = Aj[ffc] +afc_|_i(hj(V[£fc]) — Ui[£fc] +dj). Thus, © 
coincides with standard dual gradient method in (IT6b . and 
the convergence results in lf35l Prop. 8.2.6], J37|j carry over 
to this ideal setup. In this work, convergence of the system 
outputs {yi{t)}i & jy D to the solution of (PI) is assessed in 
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Fig. 2. (a) Conventional distributed optimization setup : Problem 0 is solved in a distributed fashion by using steps (H; once the problem is solved 

(i.e., iterates in {Tb} have converged to the optimal primal and dual values), the optimal reference signals {u° pt }iejv D are dispatched to the PV-inverters. 
(b) Proposed optimization-centric control architecture : The discrete-time control signal Uj [£&] generated by the dynamic controller i E Md is dynamically 
applied as an input to the inverter by utilizing a sample-and-hold (S/H) unit; the instantaneous inverter output is sampled and utilized for updating the control 
signals. The same architecture is utilized in the asynchronous case ED. upon substituting ffScil with pa . As claimed in Theorem 1, the inverter outputs 
{y;(f)}i6JV D converge to the solution of the OPF problem. Details on the distributed implementation are provided in Section Hvl 


the more general case where update of reference signals 
may be performed faster than the systems' settling times and 
asynchronously, in order to achieve the following operational 
goals: 

(01) Instead of waiting for the underlying systems to converge 
to intermediate reference levels {u i[tk]}iejy D , steps (1 1 8bb 
(fTM are performed continuously (within the limits of af¬ 
fordable computational burden); i.e., at each instant tk , one 
may have that lim^j- ||y,(f) — Ui[ifc]|| ^ 0 for at least one 
dynamical system. This scenario is particularly relevant since 
step ( 1 1 8cb is computationally light: it affords a closed-form 
solution when the inverter is operated under (cl) and (c2), 
and it involves a projection onto the inverter operating region 
under (c3) [cf. Figure [T). 

(02) The computational time required to solve the SDP prob¬ 


lem (1 1 8db is typical higher than that required by the projection 
operation (flScil : especially when (1 1 8cb affords a closed-form 
solution (see e.g., ED, and pertinent references therein). Thus, 
convergence of the system outputs is investigated for the case 
where the update of the input reference levels {u;[ffc]}ieA/D 
and the dual variables {\[tk]}iejy D is performed at a faster 
rate than ( 1 1 8db . 


To this end, suppose that the computational time required 
to update matrix V spans M < +oo time intervals; that is, if 
the computation of ( 1 1 8db starts at time tk based on the most 
up-to-date dual variables {Xi[tk]}ieAT D , its solution becomes 
available only at time tk+M- In contrast, the controller affords 
the computation of steps (1 1 8cb and (1 1 8bb at each time {tk}ke n- 
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To capture this asynchronous operation, consider the mapping 



c(k) := M 

k 

k £ N. 

(19) 

Using (IT9l>. steps (fl8l> for all i £ 

Md are modified as: 


Yi [tk] 

= r i^x i (f fe ),dij 



(20a) 


= A i [t k ] T cx/c+i 

i( v [4(fc)D - yi[tk] + d i) 

(20b) 


= arg min G, ; (u !: ) 

u i£y>i 


(20c) 


for all tk, k £ N. Further, matrix V[t c (j.)] is updated (at the 
possibly slower rate) as: 

V[i c(fc) ] = arg min H(V) T ^ Aj[f c(fe) ] h;(V). (20d) 

Since c(fc) = k over the interval {tk, ■ ■ ■, tk+M-i}, (120db 
indicates that V is being updated every M time slots. The 
block diagram for (l20l > can be readily obtained by replacing 
step (1 1 8db with (120db . as well as (1 1 8bb and (1 1 8cb with (120bb 
and ( 120cb . respectively, in Figure 0 

In the following, convergence of the system outputs to 
the solution of the steady-state optimization problem (PI) is 
established when the reference signals are produced by (l20l >. 
Of course, by setting M — 1, steps (l20l > coincide with ( ITSl ). 
and therefore the convergence claims for this more general 
setting naturally carry over to the synchronous setup in (fl8l >. 

For brevity, collect the system outputs in the vector 
y := [y|,..., y^J 1 , and the dual variables in A := 

[A{,..., A^- D | T . In the following, it will be shown that (1 1 8bb 
and (120bb are in fact e-subgradient steps 12T1 Proposition 2] 
whenever liin, t - ||y(f) — u[i/ c ]|| ^ 0 and/or M > 1. Before 

°k 

elaborating further on the error e[t k ], notice that from the 
compactness of sets V and {T’^ieAT^ it follows that there 
exists a constant 0 < G < +oo such that the following holds: 

|| h( V) — y + d|| 2 < G, VV€V,Vy € ^ (21) 

with y := )^i x )^2 x ... x y^. Furthermore, given the 
Lipschitz-continuity of the contraction mapping ( 1 1 8cb Il44l 

tq(A) = arg min Gi(uj) — A T u,, Mi £ Md (22) 

u i^yi 

there exists A [tk] satisfying 

y i[tk\ = arg min Gi(uj) - a}[4]u,, Mi £ M D (23) 

that is, yi [tk] would be obtained by minimizing the Lagrangian 
L(V, u, A[t fc ]) when A [t k ] := [AJffc],.. •, A^ffc]] 1 replaces 
A [tk]- The following will be assumed for A[ffc], 

Assumption 4: There exists a scalar G, 0 < G < Too, such 
that the following bound holds for all tk, k > 1 

IIA [tk] — \[t k ]h < G||A[ffe] — A[ffc_i]|| 2 . (24) 

Condition (l24t implicitly bounds the reference signal tracking 
error ||y[ffc] — u[ffc]|| 2 , as specified in the following lemma. 


Lemma 1: Under Assumption!!] it follows that the tracking 
error ||y[ffc] — u[ffc]|| 2 , k £ N, can be bounded as 

\\y[t k ] - u[t k ]\\ 2 < LGGa k (25) 

where L is the Lipschitz constant of function u, (A) in (l22l >. 
Proof. See Appendix lAl □ 

It can be noticed from (l25l > that the tracking error is allowed 
to be arbitrarily large, but the outputs y[t k \ should eventually 
follow the reference signal u[tfc], In fact, since the sequence 
{a k } is majorized by {rj k }, and rj k ! 0, it follows that 
l|y[*k] - u[f fc ]|| 2 —> 0 as k —> oo. Based on this assumption, 
two results that establish convergence of the overall system are 
in order: Lemma \2} provides an analytical characterization of 
the e-subgradient step, while 77;eore;n|7]cstablishes asymptotic 
convergence of the output powers to the optimal solution 
of (PI). 

Lemma 2: Suppose that at least one of the following state¬ 
ments is true: i) M > 1; ii) at time t k , yi[t k \ ^ u*[tfc] for at 
least one dynamical system. Then, h(V[f c ( fe )]) — y[t k ] + d) 
is an e-subgradient of the dual function at A [t k ]- In particular, 
under Assumption 0 and with M < +oo, it holds that 

(h(V[t c ( fc) ]) - y [t k ] + d)) T (A - A [t k ]) 

> q( A) - g(A[ffc]) - e[t k ] V A (26a) 

where the error e[t k \ > 0 can be bounded as 

k—c(k) 

e[ffe] < 2 a k GG 2 + 2G 2 ^ ] a k -h +1 • (26b) 

h=l 

Proof. See Appendix |B1 □ 

Theorem 1: Under Assumptions 00 and for any 1 < M < 
Too, the following holds for the closed-loop system (l20l > when 
a stepsize sequence {oifc}fc e N satisfying conditions (sl)-(s3) 
is utilized: 

(i) Aj[ffe] —x A° pt as k —> oo,Mi £ Nd\ 

(ii) V[t c(fc )] -x V opt and {u-x u° pt } ieA x D as k -x oo; 

(iii) y i(t) —X u° pt as t — X oo, V* £ Md- 

Statements (i)-(iii) hold for any initial conditions 
V[0], {u, ; [0]} iej v- D , {yjCO)},^^; {A,;[0]} iej y D > and any 

duration of the intervals 0 < t k — t k -i < oo, k £ N. 

Proof. See Appendix ICl □ 

Remark (convex relaxation or approximation of the OPF). 
For illustration purposes, the SDP relaxation technique for the 
OPF task is considered in this paper. However, the synthesis 
procedure outlined in the next section to develop feedback con¬ 
trollers that drive the inverter outputs to solutions of pertinent 
convex optimization problems can be utilized in a variety of 
different setups. For example, it can be utilized when second- 
order cone relaxations [251 or linear approximations II261 - H28I 
of the OPF problem are utilized. The paper considers convex 
relaxations or approximations of the OPF problem because 
dual e-subgradient-type methods are guaranteed to converge 
to optimal dual and primal solutions when applied to convex 
problems. The design of feedback controllers in the case of 
non-convex programs and their convergence will be the subject 
of future efforts. 
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Algorithm 1 Distributed architecture: inverter operation 

for k = 1, 2, 3,... do 

[51] Sample the inverter output y j[t fc ] = [Pi[tk], Qi[tk]] J ■ 

[52] Receive hi(V[t c (fc)]) from utility, if available (i.e., if 
c(k) = k). 

[53] Compute stepsize otk+ i, and update via 120bt . 

[54] Update the setpoints Ui[£fc+i] via ( I20ct . and implement 
Ui[t fc+ i] at the inverter. 

[55] Transmit Aj[£fc+i] to the utility. 

Go to step 1. 

end for 


Algorithm 2 Distributed architecture: DSO operation 

for k = M, 2 M, 3 M, ... do 

[51] Transmit hi(V[f c (fc)]) to inverter i. Repeat for all i £ A (d. 

[52] Receive Ai[f c (M + i] from inverter i. Repeat for all i £ Md. 
[S2] Start the update of V via J20db . 

end for 


Remark (discrete variables). The OPF formulation considered 
in this paper does not include the optimization of the trans¬ 
former taps at the substation as well as taps of capacitor 
banks. Rather, these quantities are considered as inputs to 
the OPF problem, and are utilized to set the voltage at the 
substation 0-El and form the (time-varying) admittance 
matrix Y in (0. This strategy ensures full interoperability 
of the proposed controllers with legacy switchgear. However, 
it is worth noticing that transformer taps can be included in 
the optimization procedure by following the relaxation method 
investigated in e.g., H33I . 

IV. Distributed implementation of the 

CONTROLLERS 

When applied to the PV-inverter regulation problem outlined 
in Section III-BI the controller ( I20bb — (120db endows each PV- 
inverter i £ No with the capability of steering its power 
output y i(t) = [Pi(t),Qi(t)\ J towards the solution u° pt = 
[P° pt , Q° pt ] T of the formulated AC OPF problem. Claims (i)- 
(iv) of Theorem [7] hold for any duration 0 < tk - tk- i < oo, 
k £ N, for any size of the distribution network. 

The feedback controller ( flol i affords a distributed imple¬ 
mentation, where optimization tasks are distributed between 
the DSO and individual PV systems; see also Figure 0 In 
particular: 

i) Updates (120bb - (120cb are implemented at each individual PV 
system (they are either embedded in the inverter microcon¬ 
troller, or, at the gateway level), and u, and A, are stored 
locally at the same inverter; these steps are performed contin¬ 
uously, within affordable computational and hardware limits. 
Particularly, (120cb is performed with the goal of pursuing 
inverter-related optimization objectives such as minimization 
of the real power curtailed 0. 

ii) at the DSO, updates (120db are performed with the goal of 
pursuing system-wide optimization objectives such as mini¬ 
mization of power losses and voltage regulation (this step is 
performed every M time steps). 

To exchange relevant control signals, a bidirectional mes¬ 
sage passing between DSO and individual PV systems is 


necessary. This entails the following message exchanges every 
M time slots: h,(V [/*-]) is sent from the DSO to inverter i; 
subsequently, the up-to-date dual variable Ai[ifc] is sent from 
inverter i to the DSO. Notice that customer i £ Md does 
not share load demand and PV-related information with the 
DSO; in fact, information about the loads is not necessary 
when computing the update ( 120db at the DSO. Exchanging 
just Lagrange multipliers rather than power iterates ensures a 
privacy-preserving operation. The operating principles at both 
inverter and DSO are tabulated as Algorithm 1 and Algo¬ 
rithm 2, respectively, and schematically illustrated in Figure [2] 
In Algorithm 1, is it also shown that the stepsize sequence 
{afc}fcgN satisfying conditions (sl)-(s3) is computed at the 
inverters’ side; when changes in the load and solar irradiation 
conditions occur (that is, the inputs of the underlying OPF 
task change), inverters exchange information to restart the 
sequence. For example, each inverter can utilize the sequence 
afc = c/s/k — n, with k > 1, c > 0 a given constant, and n 
the index of the instant t n with the last change in the operating 
conditions. 

Before proceeding, it is worth reiterating the underlying 
difference between distributed optimization approaches a, 
ED, m~M and the proposed idea. In particular: 
Conventional distributed optimization: distributed OPF ap¬ 
proaches involve the computation of steps similar to (16), 
schematically illustrated in Figure 0a). Particularly, the opti¬ 
mal reference signals u“ pt are commanded to the PV-inverters 
only after iterates (16) have converged. Accordingly, the 
reference signals are updated at a slow time scale, dictated 
by the time required to solve the distributed OPF solver. 
Proposed scheme: as illustrated in Figure 0b), the proposed 
controllers continuously update the setpoints, based on current 
system outputs, as well as solar irradiance and load conditions. 
Hence, the setpoints are updated at a significantly faster rate, 
that may be on the same order of the inverter dynamics. 
As a result, the proposed controllers dynamically refresh the 
OPF-based targets every time that there is a variation in 
loads, conventional generators, and solar irradiance, and enable 
adaptability to fast-changing conditions. 

To implement the proposed architecture, each controller at 
node i £ Md needs to collect at each time tk measurements 
of the demanded loads Pp tl and Qp., , as well as the prevailing 
solar irradiation conditions (which translate to the maximum 
available real power). On the other hand, to perform step (I20db . 
the DSO requires knowledge of the system topology and the 
load demand at nodes i £ M\Md ; of course, any AC OPF 
formulation has similar prerequisites in terms of required data 
and measurements 0-0, ED, ED, ED. Since functions 
{hi(V)}j e jvi> are li near in V, the prerequisite ( fl4l > solely 
depends on the topology of the distribution network; thus, (fl4l > 
can be checked at the utility side once matrix Y is available. 

Finally, it is worth noticing that consensus-based techniques 
can be adopted to speed up the computation of step (I20db 
and improve scalability with respect to the distribution-system 
size 0, ED, El, El- For example, by leveraging relevant 
matrix completion arguments m, Lagrangian decomposition 
and dual gradient techniques can be adopted to decompose 
the computation of (120db across lines and/or portions of the 





















IEEE TRANSACTIONS ON POWER SYSTEMS 


9 


system |f33l , |38l . The resultant algorithm would be similar to 
the one in ( l20l >. but with step (120db replaced by multiple sub¬ 
problems (one per line or portion of the system) that are solved 
in parallel, followed by relevant dual updates; see e.g., (33], 
l38l and 0, ED. It can be readily shown that the convergence 
claims of Theorem [7] carry over to this setup. 

V. Test cases 

The proposed PV-inverter control scheme is tested using a 
modified version of the IEEE 37-node test feeder and the IEEE 
123-node test feeder. The modified network is obtained by 
considering a single-phase equivalent; line impedances, shunt 
admittances, as well as active and reactive loads are adopted 
from the respective dataset^ The solver SDP3 is utilized to 
solve relevant SDPs in MATLAB, whereas the update of the 
inverter setpoints is computed in closed form. The objective 
of the test cases is to numerically corroborate the claims (i)- 
(iii) of Theorem [7] 

In the OPF problem, the voltage limits V min and l /max are 
set to 0.95pu and 1.05pu, respectively. In the first test, the 
IEEE 37-node test feeder illustrated in Figure[3]is utilized. The 
voltage magnitude at the point of common coupling is fixed to 
| Vo | = 1 pu, and it is presumed that 6 PV systems are present 
in the network and they are located at nodes 4,11, 22, 26, 29 
and 32. Following the technical approach of (30l Ch. 8] 
and ED, a first-order system is utilized to model the real and 
reactive power dynamics of each inverter. Further, inverters im¬ 
plement strategy (c3), and their regions of possible operating 
points is formed based on the inverter power ratings { S, 
and the available active powers {i” v }ieA D ■ Specifically, the 
power ratings are assumed to be 50,120,50,100,120, and 80 
kVA, whereas the following values for the available powers 
p av := [P 1 av ,..., P™ D ] J are considered in order to test the 
adaptability of the feedback controller to changing prevailing 
conditions (with time intervals normalized with respect to the 
time constant r): 

(11) p av (f) = [22,67,21,50,68,40] T kW, f/r e [1,200]; 

(12) p av (f) = [25,80,24,55,85,45] T kW, f/r e [201,400]; 

(13) p av (f) = [31,92,29,65,92,54] T kW, f/r € [401,600]; 

(14) p av (f) = [26,84,25,57,86,47] T kW, f/r e [601, 700]. 
At t = 0, the output active and reactive powers are 0 kW and 
0 kVAr, respectively. No minimum power factor constraints 
are enforced (i.e., is 9 = n/2), and P™ ln is set to 0 (4|. In 
this fiesr test, 11 (V) models the cost of power drawn from the 
substation as H(V) = (Tr($ 0 V)) 2 + 10 x Tr($ 0 V). On the 
other hand, the function Gi(Pi,Qi) is set to 

Gi(Pi, Qi) = ai (Pr - Pi) 2 + bi (P av - Pi) 

+ c iQi + di\Qi\ (27) 

in order to minimize the amount of curtailed real power, as 
well as the amount of reactive power provided. It is, however, 
worth emphasizing that various alternative cost functions can 
be accommodated in the proposed framework, and is 
utilized as a representative example. Coefficients ai,bi,Ci,di 
are set to di = 1, bi = 10, q = 0.01, di = 0.01 for 

2 Available at: ewh.ieee.org/soc/pes/dsacom/testfeeders. 



Fig. 3. IEEE 37-node test feeder considered in the test cases. 

i = 1,...,4, and a* = 1,6* = 10 ,d = 0.03,di = 0.03 for 
i = 5,6. With this setup, the SDP relaxation was first tested 
with these input data, the SDP solver identified solutions with 
rank-1 matrices V opt lf23l . (241 . 

At each inverter i £ No , the reference signal u, \tk is 
updated every t = t sec; i.e., tk — tk -i = r for all k £ N. 
This implies that a new reference signal Uj[ffc] is applied to 
each inverter faster than the output power settling time (which 
corresponds to approximately 5 t for a first-order system). On 
the other hand, matrix V[tfc] is updated every t = 2rsec; 
i.e., M = 2 in (l20t . This means that the inverter setpoints 
{u,[ffe]}ig J v’D are updated at a faster rate than matrix V[£jt]. 
The stepsize in (l20l) is set to ak = d/y/k — n, with k > 1 and 
n the index of the instant t n with the last step change. 

Figure Q] illustrates the evolution of the real and re¬ 
active powers generated by the inverters. It can be seen 
that the inverter outputs {y i[tk\ = [P,;(f), Qi(£)] T } con¬ 
verge in all the considered intervals (II)—(14), and the 
output powers at convergence coincide with the solutions 
of the OPF (PI); for example, before the step change 
at t = 200r, the active and reactive powers converged 
to the OPF solution 21.8,66.9,20.9,67.9,39.9 kW and 
39.2,85.6,40.7,77.1,31.4,39.8 kVAr. This corroborates the 
claims of Theorem [[] Figure [4jb) also provides a snapshot 
of the evolution of the output reactive power for inverter 2; 
it can be seen that a new reference level is applied after r 
seconds, before Q 2 (t) settles around the intermediate setpoint. 
It is also interesting to note that, in the considered setup, the 
steady-state reactive powers coincide with the available powers 
p av (t), and reactive compensation turns out to be the optimal 
ancillary service strategy. Similar trajectories would have been 
obtained when the loads are also time varying. Future efforts 
will explore variations of load and solar irradiance that may 
have the same temporal scale of the dynamics of (l20l >. 

In the second test case, a scenario with 
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(a) 



(b) 

Fig. 4. Convergence of CD, when the inverter-power dynamics are 
approximated as first-order systems with time constant r, for four different 
solar irradiance conditions. Plots illustrate the convergence of the real and 
reactive powers to the solutions of the formulated OPF problem. 

high PV-penetration is considered. Specifically, 
17 PV systems are assumed located at nodes 
4,11,13,16,17,20,22,23,26,28,29,30,31,32,33,34,36 
of the feeder depicted in Figure 0 and their 
AC power ratings are assumed to be s = [50, 

200,220,120, 200,120,150, 50, 280,100, 250,100,120, 200, 
110,250,150] kVA. Similar to the previous case, step changes 
in the solar irradiation (and, hence, in the available powers 
p av ) are considered in order to test the adaptability of 
the feedback controller to changing prevailing conditions. 
Specifically, the following values are tested: 

(11) p av (i) = 0.7s kW, t/r £ [1,200]; 

(12) p av (f) = 0.8 s kW, t/r £ [201,400]; 

(13) p av (f) = s kW, t/r £ [401,500]; 

(14) p av (f) = 0.6s kW, t/r £ [501,700], 

At each inverter i £ Md, the reference signal u i[tk) is 
updated every t = r sec, while V[ffc] is updated every 
t = 2rsec. The voltage magnitude at the substation is set 
to 1.02 pu, while coefficients a*, 6i,Cj,di in ( |27| > are set to 
oij = 1,5,; = 1, c, = 0.5, di = 3 for all i = 1,..., 17. All the 
other simulation parameters are similar to the previous test 
case. The SDP relaxation was tested under this setup, and the 
SDP solver identified solutions with rank-1 matrices V opt for 
all the four cases considered fiH, EH. 

With this setup, when inverters operate at unity power factor 
and set Pi = P av , the voltage magnitudes exceed the upper 
limit of 1.05 pu during the interval t/r £ [401,500] in 
10 nodes. Specifically, the voltage profile is shown with the 


yellow trajectory in Figure [3 a). 

The objective of this test case is twofold: i) demonstrate 
voltage regulation capabilities of the proposed scheme; and, 
ii) demonstrate that the convergence speed is not deteriorated 
when a higher number of PV systems are controlled. As for 
objective i), it can be clearly seen in Figures 0 a) that the 
voltages are steadily kept within the limits l /n " n and V /rnax ; 
particularly, the green trajectory in Figure 0a) shows that 
the proposed scheme favors voltage regulation even during 
peak generation conditions, while minimizing the amount of 
curtailed real power [cf. (l27b l. Figures 0c) and0d) illustrate 
the evolution of the real and reactive powers generated by 
the inverters. Comparing Figure 0 with Figures 0c) and0d), 
it can be noticed that the proposed controllers still provide 
fast adaptation capabilities to changes in the solar irradiation; 
and furthermore, convergence speed is not degraded when an 
increased number of PV systems are controlled. 

In the third test case, the IEEE 123-node test feeder 
illustrated in Figure 0 is utilized. The voltage magnitude 
at the point of common coupling is fixed to |Vo| = 1 
pu, and it is presumed that 10 PV systems are located 
at nodes 15,23,47,66,71,81,86,91,108 and 110. Inverters 
implement strategy (c3), and their AC power ratings amount 
to s = [500,450,200,300,200,200,150,150,200,350] kVA. 
Changes in the solar irradiation are considered in order to 
test the adaptability of the feedback controller for this larger 
system; the following values are tested: 

(11) p av (f) = 0.8 s kW, t/r £ [1, 200]; 

(12) p av (f) = s kW, t/r £ [201,400]; 

(13) p av (f) = 0.6s kW, t/r £ [401,500]. 

Similar to the previous test cases, the reference signals 
{ui[f/c]} are updated every t = r sec, while V[ifc] is 
updated every t = 2rsec. In the OPF problem, function 
iT(V) captures (the cost of) power losses, and it is set to 
H(V) = (Tr(LV)) 2 + 5 x Tr(LV), with matrix L formed 
as described in ITOl . The coefficients in ( 1271 are set to 
a,i = 1, bi = 10, Cj = 0.5, di = 3 for all inverters. All the other 
simulation parameters are similar to the previous test case. 
The SDP relaxation was tested, and the SDP solver identified 
solutions with rank-1 matrices V opt . 

Figures 0c) and0d) illustrate the evolution of the real and 
reactive powers generated by the 10 inverters. It can be seen 
that the inverters quickly regulate the power outputs to new 
OPF setpoints. In particular, comparing Figure 0 Figures 0 
and 0 it can be noticed that the convergence speed of the 
proposed controllers is not degraded when a larger distribution 
network is controlled. Notice that inverters are required to 
curtail real power in order to adhere to voltage limits. 

Finally, it is worth mentioning that, for the update (I20db . 
well-established complexity bounds for convex programs such 
as SDPs Il43l exist; these bounds quantify how the worst- 
case computational complexity increases with the number of 
variables (i.e., the network size). Further, in case of an SDP, 
provably convergent parallelization techniques can also be 
leveraged to speed up the computation of (120db : see, e.g., 13. 
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Normalized time t/r 



(c) 

Fig. 5. Test case with high PV-penetration. (a) Voltage profile at t = 450r. (b) Real powers provided by the inverters, (c) Reactive powers provided by the 
inverters. 



Fig. 6. IEEE 123-node test feeder considered in the third test case. Blue 
squares represent nodes at which PV systems are installed. 

VI. Concluding Remarks and Future Work 

This paper considered a distribution network featuring PV 
systems, and addressed the synthesis of feedback controllers 
that seek inverter setpoints corresponding to AC OPF so¬ 
lutions. To this end, dual e-subgradient methods and SDP 
relaxations were leveraged. Global convergence of PV-inverter 
output powers was analytically established and numerically 
corroborated. Although the focus was on PV systems, the 
framework naturally accommodates different types of inverter- 
interfaced energy resources. The development of provably 


convergenct feedback controllers that seek the solutions of 
non-convex OPF formulations will be the subject of future 
research efforts. 

Appendix 

A. Proof of Lemma Q] 

Recall first that, given the strong convexity of G,(ui), it 
turns out that function u,(A) in (f22l> is Lipschitz continuous 
(in A), with a constant denoted here as L B3). Then, recalling 
that y [tk\ = Uj(A[ffc]) and u[ffc] = u,(A[tfc]), it follows that 
the left-hand side of (l25l > can be bounded as 

||y[£fc] - W[£fe] ||2 < £||A[f fe ] - A[4]|| 2 (28a) 

< LGGa k (28b) 

where (I28bb is obtained by using the following bound (which 
originates from Assumption 4): 

||A[ffc] — A[tfc ]||2 < G||A[ffc] — A[tfe_i ]||2 (28c) 

< G||afc(h(V[£ c(fc _ 1 )]) + g(y[ffc-i],d))|| 2 (28d) 

< GGa k . (28e) 

Note that (128dt follows from the dual update in ( I20bb . 
and (I28eb follows from (Oil . 

B. Proof of Lemma \ 2] 

Recall that h(V[ffc]) — u[ffc] + d is the gradient of the dual 
function (fl2l i evaluated at A [t k ] lf35l . Let g(u,d) := d — u 
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gradient of the function q u ( A) at A [t k \, and applying the 
Cauchy-Schwartz inequality, one has that 

e u [tk\ < g T (u[4],d)(A[4] - A[f fe ]) 

+ g T (y[4],d)(A[4] - A[4]) (29g) 

< 2G \\X[t k \ - X[t k }\\ 2 < 2a k GG 2 (29h) 

where (I20bb . (OH . and < l24b were used to obtain (I29hb 
from \29g) . Next, to show ( I29db . begin with the inequality 

h T (V[t c(fe) ])(A - A[f c(fc) ]) > q v { A) - qy(A[f c(fc) ]). Adding 
h T (V[f c (fc)])(A[f c (fe)] — A[ffc]) to both sides of the inequality. 


(a) 


400 r 
350 - 
300 - 
r 250 - 
^ 200 : 
C? 150 " 
100 : 
50 - 
0 - 


i> 


200 300 

Normalized time t/r 


(b) 

Fig. 7. Convergence of (20l , for the modified IEEE 123-node test feeder. The 
inverter-power dynamics are approximated as first-order systems with time 
constant r, for four different solar irradiance conditions. Plots illustrate the 
convergence of the real and reactive powers to the solutions of the formulated 
OPF problem. 


h T (V[f c ( fc) ])(A - X[t k ]) > q v { A) - qv(A[f c ( fc )]) 

+ h T (V[f c(fe) ])(A[f c(fe) ] — \[t k ]) ■ (29i) 


Adding and subtracting the sequences {<jy (A[tfc_/i+i])}f 
and {h T (V[f c ( fc )])(A[f fc _^ + i])}^“ 2 (fe) to the right-hand-side 
of (129ib . and suitably rearranging terms, one obtains 

h T ( v [4(fc)])(A - A[t fc ]) > q v ( A) - q v (\[t k ]) - e v [tk ] (29j) 

where ey\tk] is defined as 

k—c(k) 

£v[tk] := ^2 (QvWtk-hl) - Qv(^[tk-h+ 1 ])) 

h= 1 
k—c(k ) 

- ^ h T (V[f c(fc) ])(A[f fc _ ?l ]-A[f fc _ /l+1 ]). (29k) 

h =1 

Using the definition of t he gradient, the Cauchy-Schwartz 
inequality, and d2Tl> . (I29kb can be bounded as: 


for exposition simplicity, and consider decomposing (IT2t as 
<?(A) = q v ( A) + q u ( A), with 

qy (A) := min H(V) + A T h(V), (29a) 

q u { A) := min G(u) + A T g(u, d) (29b) 

u£U 

where G(u) := ff ie jy D Then, it will be shown that 

g T (y[f fc ],d)(A - A[4]) > q u { A) - q u (X[t k }) - e u [t k \ (29c) 
h T ( v [f c (fc)]) A - A[t fc ]) > qv{ A) - q v {X[t k ]) - e v [t k ] (29d) 

with e u [t k \ < 2 ct fc GG 2 and e v [t k ] < 2G 2 J2h= C i k) a k-h+ 1 - 

To show (129cb . consider the gradient of q u ( A) evalu¬ 
ated at A [t k \, which by definition leads to the inequality 
g T (y[4],d)(A - A[f fe ]) > q u (X) - q u ( A[f fc ]) for all A; then, 
add g T (y[ffc], d)(A[ffc] — A[ifc]) on both sides to obtain 

g T (y[4],d)(A- A[f fe ]) > q u { A) - q u (X[t k ]) 

+ g T (y[G], d)(A[ffc] - A[i fc ]) (29e) 

and add and subtract q u (,X[t k }) to the right-hand-side 

g T (y[ife]>d)(A-A[4]) > q u (X) - q u {X[t k ]) (29f) 
+ q u {X[t k }) - q u (X[t k }) + g T (y[ffc], d)(A[tfc] - A[f fe ]). 

In <(29B, define e u [t k } := q u {X[t k ]) - q u (X[t k ]) + 

g T (y[£fc], d)(A[tfc] — A[ffc]). By using the definition of the 


k — c(k ) 

£v[tk\ < 2G ^ ||A[ffc_h + i] — A[ffc_fe ]||2 ■ (291) 

h=l 

Finally, upon using (I20bb and d2TTl . (1291b can be fur¬ 
ther bounded as 2G|| X[t k _ h+1 ] - X[t k _ h \\\ 2 < 

o/-»2 y^k-c(k) 

2G Lh=i ' a k - h +i. 

C. Proof of Theorem [7] 

Claims (i)-(ii). Boundedness and convergence of the dual 
iterates can be proved by leveraging the results in ll22l Theo¬ 
rem 3.4]. In particular, it suffices to show that the following 
technical requirement is satisfied in the present setup: 

+oo +oo 

Oike[tk] = ^2 a k {ev[t k ] + e„[4]) < +oo . (30a) 

k =0 k—0 

From Lemma |2] it it can be shown that 

+00 +oo +oo 

Y, a ke u [tk] < J2 2a lGG 2 < 2GG 2 Y r\\ (30b) 

k—0 k—0 k= 1 

where the second inequality in (130bb follows from the fact 
that a k < rjk for all k. Since Pk < +oo, the series 
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12k=*o a k£u[tk] is finite. As for the error ey[ifc], one has that 


+oo 



-foo k — c(k ) 


y^ oLk^v [tk] 

< 

2G 2 

y^ Qfc y: &k-h+i 

(30c) 

k= 0 



k= 0 h=l 





+00 M— 1 



< 

2G 2 

y: y ak~h+i 

(30d) 




k= 0 h= 1 





+ oo M— 1 



< 

2G 2 

E E ^-a+i 

(30e) 


k=0 h =1 


where the fact that max{/c—c(fc)} = M— 1 is utilized in (130db . 
and ( I30eb follows from ( 130dl ) since the sequence {r/ k }ken 
majorizes (afcjfceN, and it is monotonically decreasing. Since 
the series {rfi,} is square-summable, a fc£y[ffe] is finite. 

Claim (iii) From the convexity of the Lagrangian in 
the primal variables, it follows that optimal primal vari¬ 
ables can be uniquely recovered as {V opt ,u opt } = 
arg minvGV.uGW L (V, u, A opt ). 

Claim (iv) At convergence, the reference signal is constant, 
with value u° pt . Then, y i(t) —> u° pt as t —> oo by (0. 

D. Extension to multi-phase systems 

For notation and exposition simplicity. Sections [TT] and [IV] 
considered a balanced distribution network. However, the 
proposed framework can be extended to multi-phase systems 
as detailed in C0|, ED, g6) and briefly explained in the 
following. 

Define as Vij C {a,b,c} and V l C {a, 6, c} the sets of 
phases of line ( i,j ) £ £ and node i £ A f, respectively. 
Hereafter, a superscript (■)'* is utilized to assign relevant 
electrical quantities to a specific phase. For example, Vf £ C 
denotes the complex line-to-ground voltage at node i £ J\f and 
phase (!) £ Vi, whereas if £ C is the phasor representation 
of the current injected at the same node and phase; further, 
Pf and Qf denote the output real and reactive powers of a 
PV-inverter connected to phase 4> £ Vi of node i £ A I'd- Lines 
(to, n) £ £ are still modeled as ^-equivalent components l32l 
Ch. 6] and the \V m n\ x \P m n\ phase impedance and shunt 
admittance matrices are denoted as Z mn £ <c\' Pm ^ y '\' PmTi \ 
and Ymh £ £\'Pm n \x\v mn \^ respectively. Three- or single¬ 
phase transformers (if any) are modeled as series components 
with transmission parameters that depend on the connection 
type (321 Ch. 8], 0. 

A prototypical non-convex OPF formulation can be readily 
obtained by enforcing the balance constraints and the voltage 
regulation constraints on a per-node and per-phase basis 0, 
ED, E3; and, by properly augmenting the cost function to 
account for (cost of) power losses and power injections over all 
phases. To develop an SDP relaxation of the non-convex three- 
phase OPF problem, consider re-defining the vector of voltages 
v as v := [vj, v|,..., v^] T , where v* := {{Vf}^ Vi ] J is a 
\Pi\xl vector collecting the voltages on the phases of node 
i £ M. Similarly, vector i now collects the currents injected in 
all nodes and phases; that is, i := [ij, ij,..., i^] 1 . As shown 
in Section [IT] Ohm’s law and Kirchhoff’s current law can be 
captured by the linear equation i = Yv where, in this case, 
the network admittance matrix Y has dimensions JA GAt IA|X 
Y^ieAf |A|, an d its entries are computed based on the system 


topology and the lines matrices {Z^, (J -) € g as specified 

in GDI, IS) and |33) . To express voltage magnitudes and 
powers as linear functions of the outer-product matrix V := 
vv H , define matrix Y f := ef(ef) T Y per node i and phase 
iT n T n T n T IT 


zttz. 

where ef := [OL, 0 

^ denotes the canonical basis of 


and { e v . 


(Pop ■ ' •> u |P<-ip e Pi ’°|P<+ip ' • ' ’°|Piv|J 


I. Next, 

£ Vi, define the Hermitian 


per node i £ Af and phase 
matrices <&f := ±(Y f + (Yf) H ), S if := §(Yf - (Yf) H ), 
and Y';' := ef(ef) T . Then, the net real and reactive powers 
injected at node i and phase tfi can be expressed as Tr(d>f V) = 
Pf — Pf i and Tr((I/fV) = Qf — Qfrespectively, whereas 
the squared voltage magnitude at the same node and phase 
reads Tr(YfV) = |V/j 2 . 

Using these definitions, the OPF problem can be formulated 
as: 


min 

v^o,{uf Gyf} 


w) + E E °t ( u ?) 

«GA/d 


(31a) 


subject to hf (V) = uf — d f, Vi £ A I’d, 0 £ A (31b) 

hf(V) = -df, Vi £ N 0 ,<j) £ Vi (31c) 

V* in < Tr(YfV) < ^ ax , V* £ A/ - , </> £ 7A (3Id) 
rank(V) = 1 (31e) 


where hf(V) = [Tr(<frfV),Tr(*fV)] T , uf = [Pf,Qf)\ 
and df = [PffjQfj 1 [cf. 0 ]. An SDP relaxation of 
problem (Oil can be obtained by discarding the rank con¬ 
straint OTB . 

The procedure outlined in Section IIII-BI can be utilized 
to synthesize controllers for the PV-inverters that solve Oil) 
in a recursive manner. To this end, if suffices to dualize 
balance constraints (Bibb to form the (partial) Lagrangian (flTT) . 
and follow steps (l20b . In particular, the resultant distributed 
algorithm entails the following operations: 

yf[4]=rf(xf(4),df) (32a) 

•Vf[tfc+i] = Af[t fe ] + a fc +i(hf(V[f c(fe) ]) -yf[t k \ +df) 

(32b) 

u f[tk+ 1 ] = arg min Gf (uf) - (\f[t k+1 ]) J uf (32c) 

ufe yf 


V[f c(fc) ] = arg min iT(V) + E E ( A f [^(fc)]) 1 hf (V) 

6V ieMo <!>eVi 


(32d) 


where (I32ab — (132cb are performed at each PV-inverter con¬ 
nected to phase (!> of node i, and (132dl) is carried out at the 
DSO. In (132dl) . the set V is defined as 


V := {V : V A 0, V^ in < Tr(YfV) < U^ ax V* £ M 
and hf (V) = -df, V</> £ V h i £ Ao} • (33) 


It can be readily shown that the convergence claims of 
Theorem [7] carry over to the multi-phase setup. 
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